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Abstract: In this publication, uncertainties in and differences between the Mc@NL0 and ^ 
POWHEG methods for matching next-to-leading order QCD calculations with 
parton showers are discussed. Implementations of both algorithms within the 04 
event generator Sherpa and based on Catani-Seymour subtraction are employed ^ 
to assess the impact on a representative selection of observables. In the case 



o 

of Mc@NL0 a substantial simplification is achieved by using dipole subtraction 

terms to generate the first emission. A phase space restriction is employed, which 

allows to vary in a transparent way the amount of non-singular radiative cor- * 

rections that are exponentiated. Effects on various observables are investigated, | 

using the production of a Higgs boson in gluon fusion, with or without an asso- 

ciated jet, as a benchmark process. The case of iJ+jet production is presented ^-H 

for the first time in an NLO+PS matched simulation. Uncertainties due to scale ^ 

choices and non-perturbative effects are explored in the production of and • ^ 

Z bosons in association with a jet. Corresponding results are compared to data knJ 

from the Tevatron and LHC experiments. 



1 Introduction 

The central topic in the development of Monte-Carlo event generators during the past decade was the 
systematic inclusion of higher order QCD effects. One of the first approaches, the merging of leading-order 
(LO) multi-jet matrix elements of varying multiplicity with parton showers, has been pioneered and matured 
in a series of papers [T], such that by now it has become the accepted standard for simulating final states 
which include multi-jet topologies. It typically yields a very satisfactory description of experimental data, but 
due to the lack of virtual corrections it can never achieve the formal accuracy of a full next-to-leading order 
(NLO) calculation. This shortcoming is especially worrisome in the case of Standard Model Higgs-boson 
production through gluon fusion, where a large if-factor indicates substantial higher-order corrections, such 
that theoretical predictions should be made at next-to-leading order accuracy or better. 

With the development of the Mc@NL0 2J and PoWHEG ;3] techniques it became feasible to combine NLO 
accurate matrix-element calculations with parton showers. This technology is called NLO matching, in 
contrast to the LO merging described above. Ultimately, NLO matching techniques allow to carry out a 
full simulation of events, including hadronisation, hadron decays and the underlying event. While Mc@NL0 



relies on a subtraction algorithm based on the parton shower approximation to collinear divergences in 
real-emission matrix elements, the POWHEG method effectively constructs a one-emission generator, similar 
to a one-step parton shower, with evolution kernels determined by ratios of real-emission and Born matrix 
elements. In this respect the PoWHEG method is very similar to traditional matrix-element correction 
techniques [1]. 

Despite having been proposed several years ago, the Mc@NL0 and PoWHEG methods were not applied to 
processes which can become singular at Born level, such as di-jet production, until recently. This delay 
signals several complications which arise in reactions with additional light partons. Di-jet production '5] and 
Z + j production 6 have now become available through the PowhegBox [7| and, more recently, plus 
two jets production was implemented [§] using Mc@NL0 methods. 

In this publication the Mc@Nlo and PoWHEG methods will be compared in some detail, and open issues 
related to their implementation and validation will be discussed. As benchmark processes, the production 
of Higgs-, W^- and Z-bosons, alone or in association with one hard jet at Born level have been chosen. 
These processes are either signals central to the experimental program at the LHC, or they contribute as 
important backgrounds to many new physics searches Apart from that, QCD-associated W^- and Z- 
production are of great interest, to study jet production and evolution in a hadron-collider environment at 
the energy frontier |101 lllj. to improve the jet energy scale determination by the experiments, or to study 
multiple parton scattering processes. At the LHC as well as at the Tevatron, typically, good agreement is 
found when comparing respective data with merged leading-order predictions or with next-to-leading order 
perturbative QCD predictions, for instance from [T^[T5] . This level of theoretical control and experimental 
precision suggests that the processes chosen here are indeed well-suited to serve as testbed for NLO matching 
methods. 

The present paper is organised as follows: In Sec. [2] basic ideas underlying the Mc@NL0 and PoWHEG 
strategies are reviewed. Similarities and differences between the two methods are discussed and potential 
pitfalls in their implementation are indicated. The solution to some of these problems is detailed in Sec. [Sj 
where the Mc@NL0 algorithm is worked out in a formalism similar to [M] . This allows to elaborate on various 
aspects of the method regarding non-trivial colour and flavour configurations. Readers not interested in the 
technical details may skip this section. In Sec. |4] perturbative uncertainties related to scale variations and 
to the different exponentiation in the Mc@NL0 and PoWHEG methods are highlighted. The impact of 
non-perturbative effects is investigated in Sec. [5] and results for W^- and Z-production are compared to 
experimental data. Finally, Sec. [6] summarises the findings of this publication and outlines some possible 
extensions of the methods presented here. 



2 Event generation with PoWHEG and Mc@Nlo in a nutshell 

This section is meant to highlight in a simple language the ideas of how NLO QCD calculations can be 
combined with subsequent parton showers. A common notation is established for PoWHEG and Mc@NL0 
and critical aspects related to the practical implementation for processes with non-trivial colour structure 
are discussed. 



2.1 Anatomy of NLO calculations 



In order to see how the existing matching algorithms work, consider first the structure of an NLO calculation. 
The total cross section can be written as 



a 



(NLO) 



d$; 



B($B)-l-V($B)-|-l(Sn*B) 



d$i 



R($fl)-D(S)($^) 



(2.1) 



where d$B and d$fl denote Born and real-emission phase space, respectively and B, V, R and D^^^ are 
the matrix elements for the Born, virtual, integrated subtraction, real emission and real subtraction contri- 
bution. Note that for hadronic initial states V is defined such as to include the collinear mass-factorisation 
counterterms. All integrands include parton luminosity and flux factors. For simplicity it is assumed that 
the processes under consideration here have no identical QCD partons in the final state, a detailed discus- 
sion of symmetry factors is postponed to Sec. [Sj In the subtraction formalisms most commonly employed 
to date [151 116j , subtraction terms can be written as the convolution of Born matrix elements with suitably 
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defined operators K, such that schematically 
d$fl D('^)($fl,) ^ d$Bd$i [b($s) ® K($i)j and I^^^^s) = / d$i B($b) ® K($i) 



(2.2) 



indicating that the two related contributions cancel locally in the Born phase space. This also suggests that 
the definition of the subtraction kernels K is somewhat arbitrary, as long as they exhibit the same singularity 
structure as the full real-emission matrix elements in the soft and coUinear limits. 

To combine an NLO calculation with a parton shower, it is essential to ensure that the result inherits the 
total cross section from the fixed-order calculation. This could trivially be achieved by multiplying a leading 
order plus parton shower simulation with a local _?C-factor - the ratio of NLO and LO cross section. The 
problem is that parton showers generate emissions which do not necessarily follow the pattern obtained 
from the real emission matrix elements, leading to a mismatch of radiation patterns at first order in the 
strong coupling constant a^. In principle, this problem has been solved by the traditional matrix-element 
correction procedure outlined in [4 , which is implemented for a number of processes in both Herwig and 
Pythia. The essence of this method lies in the fact that, for the processes it is applied to, the product of 
Born-level contribution and parton shower evolution kernel K is larger than the corresponding real emission 
term R in the complete phase space of the extra emission. This allows to correct the first (hardest) emission 
with a factor R/(BK), leading to the desired distribution in phase space. From this it can be seen that such 
an algorithm will fail for cases, where R is not always smaller than the Born-times-parton-shower result, or 
where the parton shower is not capable of filling the full phase space. In such cases a more detailed analysis 
is necessary. 



2.2 Sudakov form factors 

In this section the construction of the parton shower using Sudakov form factors will be briefly recalled. 
They are defined as 



A{t,t') = exp. 



dt 



I 77^ Jit, z, 

ZTT 



27r 



K(i, z, 



exp ■ 



d$i ^ K(<i>i) , (2.3) 



where the sequence of emissions is ordered by a parameter t - typically related to the virtuality, transverse 
momentum or opening angle of the emission. Furthermore, z is the variable defining the splitting of energy or 
light-cone momentum of the two daughters emerging in the decay of the parton, and (j) denotes the azimuthal 
angle. Collectively they define the one-emission phase space $i and its Jacobian J, while K is the splitting 
kernel encoding the detailed kinematical distribution due to soft-coUinear enhancement and spin effects. For 
simplicity, parton luminosity factors are included into this kernel, and a detailed discussion is postponed to 
Sec.[3j Note that in all practical implementations of Eq. (2.3) the argument of the strong coupling is related 
to the transverse momentum of the splitting, given by the decay kinematics, i.e. by $1. For the sake of 
clarity this argument will be suppressed in this section. 

This Sudakov form factor can be interpreted as a no-branching probability between the two scales t and t' . 
Obviously, 1 — A(t, t') then yields the probability that a splitting has taken place in the interval between the 
two scales. Therefore, to first order in as, the cross section in the parton-shower approximation reads 



4s°^ = / d^sBi^s) 



A(<o,/4) + 



d<i>i-^K($i)A(i,4) 

ZTT 



(2.4) 



where t is determined by the kinematics of the first emission, t — t(^i), such that it is smaller than the 
suitably chosen, process-dependent factorisation scale fijp, and where the arguments of the splitting kernel 
depend on this extra emission kinematics. The first term in the square bracket represents the probability 
of no extra emission to happen above the infrared cut-off to of the parton shower, while the second term 
generates one emission above this scale. It can also be seen that the square bracket as a whole integrates to 
one, exposing the probabilistic nature of the parton shower, which leaves the total cross section of an event 
sample unchanged. 
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2.3 Matrix element corrections and POWHEG 



Applying matrix element corrections to the parton shower transforms the equation above, Eq. (2.4 1, into 



JLO) _ 

'mec ~ 



d$BB($B) 



Jin £>(<1>b) 



where the modified Sudakov form factor A reads 



A{t,t') = exp- 



t' 



R($i3,$l) 



(2.5) 



(2.6) 



Again, the square bracket in Eq. (2.5) integrates to one, and the cross section of an event sample is identical 



to the Born cross section. However, the distribution of the first (hardest) emission, to first order in 
will follow the pattern given by the full real emission matrix element. It is thus correct to 0(as) if the 
upper integration limit in the real-emission term is extended to the hadronic centre-of-mass energy. This 
corresponds to a power-shower approach [T7] and ensures coverage of the full phase space. However, it also 
implies extending the resummation built into the parton shower beyond the region of its validity and is 
therefore unjustified. A more detailed discussion of this point can be found in Appendix [A] The problem 
can be solved for example in the Mc@NL0 method, see Sec. |2.4| 

In order to achieve full 0(Q;5)-accuracy in both the cross section of the produced event sample and the 
pattern of the first emission it is mandatory to replace the differential Born cross section d^^B with an 



expression that integrates to the full NLO cross section of Eq. (2.1|. This is achieved by the substitution 



B($i 



B(<I'b)+V($b) + I(S)($s)+ / d$i [r($b,$i)-D(S)($b,$i 



(2.7) 



It is straightforward to interpret this term as an NLO-weighted Born-level cross section, or, put in a slightly 

different way, a Born-level cross section that is augmented with a local /T-factor. 

Equipped with this definition the expression for the cross section in the PoWHEG-formalism reads 



(NLO) 
^POWHEG 



d$BB($B) 



A(to) 



(2.8) 



This expression is accurate up to the first order in the strong coupling constant for both the total cross section 
of the event sample and the differential cross section with respect to the kinematics of the first emission. 



2.4 Modified subtraction and Mc@Nlo 

As indicated above, hard real-emission configurations should not be included in the exponentiation in 
Eq. (2.8 1. There are also contributions which can generically not be obtained from a Born-level configuration 
by simply emitting an extra parton. Examples for this include flavour processes such as uu — cse~v^ or 
kinematical configurations such as radiation zeroes [18] . Therefore, it is well motivated to split the real- 
emission matrix elements R into an infrared-singular (soft) and an infrared-regular (hard) part, D^"^^ and 
H^'^^ respectively: R = D'"^^ -|- H^^^ |l9j. Equation (2.1) then transforms into 



fNLO) 



d$, 



B($b)+V(*b)+I^^^($s) 



d$, 



D(A)(ci>H)-D(S)($^) 



and the cross section, Eq. (|2.8|), can be rewritten as 

(NLO) 



'MC@NLO 



= / d$BB(^)($s) 



A(^)(^o) 



d$i 



d$^H(A)($^) , 
D(A)($5,$,) 



(2.9) 



B($b) 



d$«H(A)($«) , (2.10) 



where 



B(A)($b) = B($B)+V($B)+I^^H*i3)+ j d$i 



D(A)($b,$i)-D(S)($s,$i) 



(2.11) 
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and, with obvious notation, 

In this manner, the non-singular contributions at real-emission level are correctly captured and their expo- 
nentiation can be avoided. This is how the Mc@NL0 formalism works [5]. 

2.5 Subtleties in practical implementations of POWHEG and Mc@Nlo 

At the level of detail of the discussion up to now, an implementation of both methods, PoWHEG and Mc@NL0, 
seems quite straightforward. However, there are a number of subtleties, which have not been presented yet 
in a coherent fashion in a journal publication. 

1. Subtraction of infrared divergent sub-leading colour terms in Mc@NL0 

The parton shower is based on a leading logarithmic and leading colour approximation. Using it, with- 
out modification, to define modified subtraction terms D^^' will therefore necessarily miss divergences 



in sub-leading colour configurations and will result in infinite integration results in Eq. (2.11). This 
presents a clear obstacle to apply the method to processes where infrared divergences in such sub- 
leading colour configurations emerge. In [20], and referring to ideas outlined in the original Mc@NL0 
publication [2] , this problem was overcome for the case of heavy flavour hadro-production processes of 
the type gg — >■ QQ and qq — )■ QQ at Born-level. An additional factor was introduced, which multiplies 



the complete integrand in Eq. (2.11 ) and tends to zero as emissions become soft. This same method is 
applied in the case of the recently presented aMc@NLO [HI [52] . 

We propose a different, exact and process-independent solution here, which relies on choosing D^^^ 



D^^-*, leading to a tremendous simplification of Eq. (2.11 ). In the remainder of this publication we will 



refer to this scheme as the McQNlo method. Technical details of how D^^^ is exponentiated into a 
Sudakov form factor will be given in Sec. [3j 

2. Choice of scales in PoWHEG and MC@Nlo 

It was argued in [3], based on [23], that the proper choice of scale for the strong coupling in the 
kernel R/B of the PoWHEG Sudakov form factor is given by the transverse momentum of the emission 
generated in the branching process. Similarly, it was pointed out in [5] that the scale in the evolution 
kernel of an Mc@NL0 should be the transverse momentum. This particular choice of scale effectively 
resums a certain class of universal higher-order corrections. It typically leads to a very good agreement 
between the results from existing parton shower algorithms and experimental data and is therefore 
employed in all standard Monte-Carlo event generators [24j . 

Such a scale choice ultimately implies higher-order corrections to the event kinematics. The key point 



is that the differential NLO cross section defined by Eq. (2.7) can be evaluated with arbitrary scales. 
However, the functional form of the scale must be infrared and coUinear safe. This is not the case 
for the transverse momentum in the branching process, because, in contrast to the Sudakov branching 
algorithm, there is no infrared cutoff in the calculation of B. Thus, the first-order expansion of the 
Sudakov form factor does not reproduce the exact same real-correction / subtraction terms as the 
NLO calculation. The differences are numerically large in the logarithmically enhanced regions of 
phase space, where it is not expected that PoWHEG or Mc@NL0 reproduce the exact leading-order 
result. But even the region of hard radiation can be affected, depending on how much phase space is 
covered by D*^^^. 

3. Exponentiation of non-leading logarithms in PoWHEG 

Ultimately, there is a last subtlety, which has also been discussed in [7 , summarising some previous 
work by the same authors. The evolution kernel R/B in the PoWHEG method generates subleading 
logarithms and it is somewhat questionable in how far they should be exponentiated. The difference 
compared to evolution kernels in the parton shower may lead to sizable effects, up to orders of mag- 
nitude, when results for the hard, non-logarithmic tails of distributions in the PoWHEG approach are 
compared with those of pure NLO calculations, see for instance [53]. 

In this context, it is important to distinguish two effects: One is exponentiation of R/B beyond the 



factorization scale fip (cf. Eq. (2.8)), the other is the difference in the functional form of the evolution 
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kernels itself. The former effect can be emulated in our implementation of the MC@Nlo formalism and 



we will examine it more closely in Sees. |4.1| and 4.2 We use CS subtraction and a parton shower built on 
CS subtraction kernels |26] . Initially the starting scales of the parton shower will be maximised for the 
first emission, in order to fill the full phase space, similar to the POWHEG method and a power-shower 
approach [T7]. To show that the main difference between PoWHEG and Mc@NL0 lies in the enlarged 
phase space in POWHEG, a simple phase-space constraint is added to the subtraction |57j , limiting the 
logarithmically enhanced region of the real-emission contribution and adding the hard regions to the 
regular contribution. The impact of the choice of this constraint is analysed quantitatively in Sec. [4] 
The respective variations in the result are referred to as "exponentiation uncertainties" . It is important 
to note that this uncertainty is minimised in the Mc@NL0 method, as the phase space constraint is 
determined by the resummation scale. In this publication we simply use the Mc@Nlo framework to 
make the problem in PoWHEG explicit. 

3 PowHEG and Mc@Nlo in detail 

When combining next-to-leading order matrix elements with parton showers, the logarithmic corrections 
encoded in the Sudakov form factor must be matched to the full next-to-leading order prediction of the 
parton-level calculation. This is usually achieved either by subtracting the first-order expansion of the 
resummed result from the NLO calculation (Mc@Nlo with D^^^ = D*^^'), or by exponentiating the full 
real-emission correction in the resummation (Powheg). One can also construct mixed schemes, as argued in 



Sec. 2.4 The difference which is induced by the exponentiation of subleading corrections vanishes for soft or 
coUinear parton emission. However, it might play an important role in other regions of the phase space. In 
the following, a formalism that allows to discuss the difference between the Mc@NL0 and PoWHEG methods 
in more detail is reviewed. 

3.1 Notation and definitions 

Denoting sets of n particles in a 2 — )• (n — 2) process by {a} — {ai, . . . , a„}, while their respective flavours 
and momenta are specified separately through {f} — {fi,---,fn} and {p} = {pi,...,p„}, the generic 
expression for a fully differential Born-level cross section can be written as a sum over all contributing 
flavour configurations: 

dfTsCM) - ^daB({S}) , where dfTB({a}) = d$B({p}) B({a}) , (3.1) 
{/} 

The individual terms in the sum are given by 

B({a}) = C{{a})B{{d}) , 6({a}) - jt}^^ -4- IMsf {{S}) , 

\\Pf ) i>({j I) ^2 2^ 

dci>B(M) = ^^d^>B(M) , C{{d};f^^) =Xi//,(xi,m') X,ff,{x2,f^') , 
X'^ X2 

where ({q-}) denotes the partonic matrix element squared, d<?s({p}) is the corresponding differential 

n-particle partonic phase-space element, S{{f}) is the symmetry factor, F{{p}) is the flux factor, and C is 
the parton luminosity 

In a similar fashion, the real-emission part of the QCD next-to-leading order cross section can be written 
as a sum, depending on parton configurations {ai, . . . , A„+i}, by replacing the Born-level matrix elements 
B with the real-emission matrix elements TZ and the Born-level phase space d^g with the real-emission 
phase-space d^ji. 

It is then useful to introduce a notation for mappings from real-emission parton configurations to Born-level 
parton configurations and vice versa. They are given by (cf. |14j ) 



^ In the case of leptonic initial states, and ignoring QED initial-state radiation, the parton distribution functions /(a;,^^) 
are replaced by 5(1 — x). 
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The map bij,k{{A}) combines the partons Ai and Aj into a common "mother" parton a^j, in the presence of 
the spectator Ak by defining a new flavour /jj and by redefining the particle momenta. The inverse map, 
r~ J. ({a}) determines the parton configuration of a real-emission subprocess from a Born parton configuration 

and a related branching process fc — > ij, k. The radiative variables ^'^j^, denoted $i in Sec. j2jfor brevity, 
are thereby employed to turn the n-parton momentum configuration into an (n + l)-parton momentum 
configuration. 

The real-emission matrix elements, TZ{{a}), can be approximated in the soft and collinear limits by subtrac- 
tion terms 2?^^ ^.({a}), which capture the universal singularity structure when two partons i and j become 
collinear, or one of them becomes soft in the presence of a spectator parton k. 

ij collincar 

^ .^ft p(S)({^}). (3.4) 

These subtraction terms are related to the ones defined in Sec. ^through D(s) ^ Ci{A})vl^l{{A}), with 
implicit notation of dipole indices, phase space and flavour configuration on the left hand side. They are not 
uniquely defined, but can be constructed, for example, using the Catani-Seymour method jl5j or the FKS 
approach |16j . Furthermore they can be restricted in phase space 27 or extended with an arbitrary finite 
function. Similar terms can also be computed in a parton-shower approximation as 

n'^^^-V^!>{{A}) = Bib^.Mim ^^'''f,%^^^ -^8na.lC.,MP^,P^,P^) , (3-5) 

where lCij,k denote the parton-shower evolution kernels. The parton-shower approximation is meaningful 
only in the collinear region, as it implements an exact factorisation of the colour structure into a Born part 
and an emission part. In many cases also a factorisation of the helicity structure is assumed. 



3.2 From fixed-order to resummation 



Section [2] only introduced expressions for the total cross section in the MC@Nlo and POWHEG methods. A 
detailed discussion requires an analysis of the expectation value of arbitrary infra-red safe observables. The 
respective formulae are developed in the following. 

Using the subtraction terms introduced in Eq. (3.4), the expectation value of an observable O is determined 
to next-to-leading order accuracy as 



(NLO) 



= J2j d<i>s(M) 
{/} 



B({a})+V({a})+ 



o(M) 



E / d<i>fl({p}) 



Ri{A})Oi{p})- Dlf,({A})0(6,,,fe({p})) 



(3.6) 



where l'-^-'({a}) represent the subtraction terms D^^-'d^}) integrated over the extra-emission phase space. 
Note that the configurations {f}, {p} and {a} on the second line each include one more particle. 

This equation can be modified by adding and subtracting an additional arbitrary set of subtraction terms 
D*^^-'({a}) with the same kinematics mapping as D(^^({a}) 



(NLO) 



d<i>s(M)B(^)(M)0(M) 

{/} 



E / d<i>«({p}) 

{F} 



Ri{A})Oi{p})- E DLi({^})0(6.,,.({p})) 



(3.7) 



where the function B^^^da}) is defined as 

B(A)({a}) =B({a})+V({a})+ ^ Ir'i(M) 



{»l,fc} 



E E 

{Ij,k} }i=1,9 



(3.8) 
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and where d<I>^j^ represents an integral over the radiative phase space. 

When combining fixed-order calculations with resummation, the task is to define a unique starting condition 
for the parton shower. As was argued in [2], it is not possible to process the terms R({a}) and D^f 



in Eq. (3.7) separately because this would lead to double-counting. Instead, the problem is solved if the 



observables on the second line of Eq. ( |3.7[ ) are both assumed to depend on the momentum configuration 
{?}. This leads to 



(O) 



(NLO) 



J2 j d<i>Bm)^^^\{p})om) 
{/} 



J2 I d*fl({^}) 

{F} 



Ri{A})- J2 Dli({^}) 



(3.9) 



0({p}) + (0)(^°") , 



introducing the correction tcrm[^ 



{F} 



o{{P})^o{h,A{F})) 



(3.10) 



The essence of both the Mc@NL0 and POWHEG methods is to generate Eq. (3.10) by using a Sudakov 



branching algorithm, which is formulated in terms of an evolution variable t, where t cc 2piPj. Let a 
corresponding form factor be defined as 



(3.11) 



where 



A^it; {a}) = cxp <( - 5: / d<i>;^j^ e(t($^|^) - t) 



F,=q,g 



1 S{T,^;,{F,-{f})) Dgi(r~,s(^„$|^;{a})) 



(3.12) 



S,. 



s{{n) 



ma}) 



The ratio of symmetry factors, including Sij, is explained in detail in |14| . It accounts for the way in which 
the phase space is successively filled by a parton shower. 

It can then be proven fl\ (31, [TIJ , that the following formula for the expectation value of an infrared safe 



observable reproduces Eq. (3.9) to 0{as): 



(NLOMC) 



= ^1 d$B(M)B(A)({a}) 
{/} 



A(A)(to;{a}) Om) 



unresolved 



E E /d*^ii^0W*;^iB)-io)O(r,,,,(M)) 



1 DS(^i,fe(^«><i^KiB;{«})) 



S{{f}) 



B({a}) 



A(^\t;{d}) 



resolved, singular 
(A) 



i'j,k 



o({p}) 



(3.13) 



resolved, non-singular 



NLO Monte-Carlo events reproducing Eq. (3.13) can be generated in the following way: 



^ The dependence of the observable O on the kinematics of the partonic final state makes the need for a correction term 
^Q^(corr) manifest. This term integrates to zero in Eq. l|2.9[|, as 0=l and the phase-space dependence is trivial. 
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• A seed event is produced according to either the first or the second line of Eq. (3.9), 
not including the correction factor (O) . 

• If the second line is chosen (H-event) the event has real-emission kinematics and is kept as-is. 
This generates the "resolved, non-singular" term of Eq. (3.131 



• If the first line is chosen (S-event), the event has Born-like kinematics and is processed through a 
one-step Sudakov branching algorithm described by Eqs. (3.111 and (3.12). The necessary techniques 



for that case will be detailed in Sec. |3.4[ An emission might or might not occur, which is represented 



by the "resolved, singular" and "unresolved" terms of Eq. (3.13), respectively. 



3.3 The Mc@Nlo and POWHEG methods 



Equation (3.13) becomes particularly simple when 



(3.14) 



In this special case the integral over the radiative phase space in Eq. (3.8) is zero. We will refer to this as 



the MC@Nlo technique. It can be thought of as a method which uses the Sudakov branching probability for 
subtraction or, conversely, which uses the subtraction kernels for parton evolution. Possible event generation 
techniques to achieve this will be discussed in Sec. |3.4| 

The PoWHEG method, in its original form, on the other hand, defines 



^'ifl^Pr3M{A})Rm) where 



(3.15) 



It can be thought of as exponentiating the full radiative corrections into a Sudakov form factor. As such, it 
bears strong similarity to matrix-element corrected power showers. This aspect has been discussed extensively 
in [H El- 



One can construct a mixed scheme, where T>\^], is defined as 

[r({a})-rW(W)' 



D 



(3.16) 



with R*^''-' an arbitrary infrared-finite part of the real-emission cross section and Pij^fc given by Eq. (3.151 



This method leads to the original Mc@NL0 prescription, proposed in [2|, if = Dif^"*- It was also used 



,(PS) 



ij.k 



in [T9j to deal with the problem of radiation zeros in the W^* -production process in the PoWHEG approach. 

f A) 

Note that all the above choices of ^ fulfil the requirement that the kinematics and fiavour mapping be 



identical to the one in D^^l, cf. Eq. (13. 71) 



3.4 Event-generation techniques 



In the parton-shower approximation, Eq. (3.12) would read 



Ar,l(i',^";M)-cxp - 



Fi=q,g 



dt 

T 



1 a. 



27r 



Jij^k(t,Z,(j)) 



bij Zn 



(3.17) 



where z is called the splitting variable of the parton-shower model and Jij_k is the Jacobian factor associated 
with the transformation of d$^j^ into dtdzdcj) |14j . 



It is well known how to generate emissions according to Eq. (3.17). The task of a proper implementation 



of Mc@NL0 and PoWHEG is, however, to employ Eq. (3.12). In our formulation of Mc@NL0, i.e. with 



D 



(A) 



ij,k 



(S) 

^ij this can involve the exponentiation of subtraction terms which are negative, due to subleading 
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colour configurations in the real-emission matrix elements, cf. Sec. |2.5[ This leads to a form factor larger 
than one, which cannot be interpreted in terms of a no-branching probability. In this section we will show 
how to generate such form factors efhciently in a Markov chain Monte Carlo. The respective algorithm is 
the basis for our implementation of the Mc@NL0 method. It was outlined in [25] and reformulated in [3D]. 
As it needs to be modified slightly in the current context, it is briefly discussed here again. The method is 
based on an extension of the veto algorithm ^31j . 



The veto algorithm 



Let t be the parton-shower evolution variable and f{t) the splitting kernel JC, integrated over the splitting 
variable The differential probability for generating a branching at scale t, when starting from an upper 
evolution scale t' is then given by 



V{t,t') = fit) exp. 



dtfit) 



A new scale t is therefore found as 
t ^F-'[F{t')+\og4] 



where 



dtfit) , 



(3.18) 



(3.19) 



and where # is a random number between zero and one. The key point of the veto algorithm is, that even if 
the integral Fit) is unknown, one can still generate events according to V using an overestimate git) > fit) 
with a known integral Git). Firstly, a value t is generated as t = G^^ [Git') + log#]. Secondly, the value 
is accepted with probability fit) /git). A splitting at t with n intermediate rejections is then produced with 
probability 



Vnit,t'\ 



m 

9{t) 



git) exp 



n 



digit) 



dt,: 



g{tt) 



exp 



digit) 



(3.20) 



where in+i = t' and to = t. The nested integrals in Eq. (3.20) can be disentangled, and summing over n 



leads to the exponentiation of the factor git) — fit), such that Eq. (3.181 is reproduced [3T] . 



Analytic weights 



The purpose here is to introduce an additional overestimate /i(t). The additional weight git) /hit) is then 
applied analytically rather than using a hit-or-miss method. This leads to the following expression for the 
differential probability to generate an emission at t with n rejections between t and t' 

dt/i(t)| 

(3.21) 



rnit,t') =^hit) exp 
9{t) 



n 



dt,, 1 



m 

9{ti) 



tiiti) exp 



dthit) 



Note that the function hit) and the ratio fit) /git) must be positive definite. Generating events according to 



Eq. (3.21 1 is then straightforward. In order to recover the desired distribution, the following analytic weight 



needs to be applied 

wit,ti, ...,tn) = 



9{t) -Q 9{t 



hit^) - fit,) 



hit) f \ hit,) git,) - fit,) 



(3.22) 



Here the term git) /hit) is due to the acceptance of the emission. The product, which is needed for an 
exponentiation of hit) — fit) instead of git) — fit), runs over all correction weights for rejected steps. 



Equation (3.22 1 can lead to negative weights, which reflect the fact that sub-leading colour configurations 



are taken into account and that the a-priori density hit) might underestimate fit). 

For simplicity, the existence of only one splitting function is assumed, i.e. that there is no flavour change of the splitter 
during the evolution. The extension to flavour changing splittings is straightforward, but it would unnecessarily complicate the 
notation at this point. 
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Implementation of the MC@NlO method 

(A) (S) 

In order to implement the Mc@NL0 method with D^^ = D^-^- ^ we need to identify the function / above 

with the (z, 0)-integrated subtraction terms f., divided by the Born contribution. A convenient choice of 

the function h will be the (z, 0)-integral of the parton-shower evolution kernels Dj-J^\ divided by the Born 
contribution. We are then free to choose the auxiliary function 5 on a point-by-point basis, but a convenient 
way is to define g = C f, where C is a constant larger than one. This guarantees that both acceptance and 
rejection term are generated in sufficient abundance to reduce Monte-Carlo errors. In the remainder of this 
paper, the functions /, g and h are thus schematically identified as follows: 

/^D(^), /i^D(ps), and g^Cf, (3.23) 

where C « 2 is a constant. 

is) 

By exponentiating Dj^- j, in this manner, we can suppress the integral over the radiative phase space in 



Eq. (3.8) and thus substantially simplify the event generation. 



4 Analysis of perturbative uncertainties 



In the following the uncertainties discussed in Sec. 2.5 are investigated in detail. The event generator 



Sherpa 02] sets the framework for this study, including its automated Mc@NL0 implementation, the 
matrix-element generator AMEGIC+-1- [33] , an automated implementation [34] of the Catani-Seymour dipole 
subtraction method [IF and the parton shower model described in [^SlUn]- The CTEQ6.6 PDF set 35J is 
used together with the corresponding parametrisation of the running coupling. All analyses are carried out 
with the help of Rivet [35] . 

4.1 Higgs-boson production in gluon fusion 

The production of a Higgs boson at the LHC with 7 TeV centre-of-mass energy serves as a first example. This 
process amplifies some of the effects discussed in Sec. |2.5| because the radiating partons are gluons. There 
are thus no valence PDF's involved at Born level, which allows to test fixed-order effects and resummation 
in a relatively clean setting. Results presented in this section do not include non-perturbative or multiple 
parton scattering effects. A detailed study of the POWHEG algorithm and its associated uncertainties has 
already been presented in fMl. In the following, the focus will therefore mainly lie on our implementation 
of the Mc@NL0 method, but comparisons with PoWHEG will be made at the appropriate junctures. 

For an unbiased comparison of the Mc@Nlo and PoWHEG techniques, the events are analysed with a 
minimal set of cuts. A Higgs-boson mass of nih — 120 GeV is assumed and two r-leptons with |?7| < 3.5 
and p± > 25 GeV are required; jets are defined according to the inclusive fcj^-algorithm [37] with R = 0.7 
and p_L.min — 20 GeV. The renormalisation and factorisation scales are set to fip = fJ-R = mh- The effective 
coupling of the Higgs to gluons, mediated by a top-quark loop, is modeled through an effective Lagrangian 
|38j . Both powers of as in the effective gluon-gluon-Higgs coupling are also evaluated at m/j. 

The comparison with fixed-order predictions and predictions from standard parton showers (referred to as 
L0(8)PS in the following) presents a first crucial test of the MC@Nlo and PoWHEG methods. In this context, 
Sherpa is also used as a framework for NLO fixed-order event generation, enabling a comparison of all 
approaches with identical input parameters. In order to vary the amount of exponentiated real-emission 



corrections in our MC@Nlo, which is governed by Eqs. (2.10)-(2.12) or, more rigorous, Eqs. (3.12|-(3.14|, 
the phase-space restriction described in [37] is employed and the parameter acut is varied. It has the effect 
of increasing (acut — > 1) or decreasing (acut — > 0) the phase space for non-singular contributions in D'^^^ 
by setting an upper bound on the virtuality t of the splitting parton. The leading logarithm is then of 
the form as log^(t/acut s), with s the centre-of-mass energy of the colliding protons. Hence, for acut ~ 1 
the leading logarithm contains the wrong argument acuts ^ ^J■'p, clearly at odds with usual choices of 
the resummation scale and in violation of factorisation theorems. Similarly, the original PoWHEG method 
exponentiates radiative corrections throughout the real-emission phase space, corresponding to our MC@Nlo 
with acut — 1- Some practical implementations of PoWHEG propose to suppress non-singular terms with a 
continuous function which tends to one in the singular limits and approaches zero in those regions of phase 
space where emissions become hard [TH] ■ However, even in this approach the volume of phase space where 
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Figure 1: Transverse momentum of the Higgs boson in inclusive Higgs boson production (m/j — 120 GeV) 
at Ecms — 7 TeV. The variation of MC@NLO predictions with varying acut (denoted a in 
the legend) is shown in Fig. [(a)] while Fig. [(b)] compares the MC@Nlo, Powheg and LO(g)PS 
methods. 



exponentiation is performed is left invariant, and only the influence of non-logarithmic terms is reduced. The 
leading logarithm is thus still of the form as log^(i/s). The traditional LO(X)PS method directly implements 
the DGLAP resummation and thus uses the factorisation scale fip as a phase space constraint, resulting in 
the correct logarithmic form log^(i/^|n). It is clearly desirable to implement this very same constraint 
in our Mc@NL0. This requires the computation of a new class of integrated subtraction terms in the 
Catani-Seymour framework, which is beyond the scope of this work and will be the topic of a forthcoming 
publication. 

The influence of ctcut on the MC@Nlo predictions and its relation to the PoWHEG and LO(X)PS predictions 



are investigated for various observables in Figs. T][5 Taking the 0{as) expansion of Eq. (3.13) at face 



value leads to the expectation that results from Mc@NL0 (for any Ofcut) and PoWHEG should approximately 
agree with flxed-order predictions in the limit of hard, well separated partons. In the limit of soft/collinear 
radiation the common logarithmic structure of the resummation in Mc@NL0, Powheg and LOiSjPS should 
lead to identical results in all three approaches. 

In practice, the transverse momentum spectrum of the Higgs boson, shown in Fig. |la[ exhibits a large 
dependence on acut- This roots in the different composition of the distribution in terms of H- and S-events, 
which is depicted in Figs. [6] and [7] Neglecting subsequent parton shower emissions, H-events above acut 
coincide with those of the NLO calculation, while S-events resum the leading logarithm of the real emission 
matrix element and occur with an additional enhancement B'-^^'/B. This enhancement, though formally 
contributing to C(a^), can be numerically large. Thus, when allowing S-type events to flU more than the 
phase space suited for the resummation of the DGLAP logarithms, their emission probability is not only 
distorted by the exponentiation but also by an artiflcial if-factor B^^^/B, as compared to the NLO result. 
For smaller values of acut, the argument of the leading logarithm is closer to the one in DGLAP evolution, 
even though acut has a functional form which makes a direct comparison difficult. It is thus observed that 
the MC@Nlo distributions come close to the flxed-order result in the hard region, and the method behaves 
exactly as expected. 

Further, excessively small values of acut lead to spurious results, as the S-event prefactor B'^^^/B turns 
negative. This is exemplified in Fig. |6j Thus, those values of acut are preferred which ensure that S-events 
only fill the phase space suited for resummation and simultaneously avoid 

B(A) 

turning negative. Defining 

such an allowed range of acut introduces an unwanted process dependence into this implementation. Clearly 
a physically more meaningful definition of phase space constraints will ultimately lead to removing this 
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Figure 2: Transverse momentum of the Higgs boson (left) in inclusive Higgs boson production (m/i = 120 
GeV) at Ecms = 7 TeV. The soHd Hues show the full Mc@NL0 (green) and POWHEG (red) results, 
while the corresponding dashed lines show the results where the prefactor of the resummed S 
events, B and B^^^ respectively, have been replaced by B. In all cases only the first emission is 
taken into account. The black dashed line shows the NLO fixed-order prediction. 



uncertainty. 

Comparing these results to the predictions of the PoWHEG and LO(g)PS methods, depicted in Fig. lb one 
finds that POWHEG behaves as our Mc@NL0 implementation with acut = 1, suffering from a too large phase 
space for exponentiation. On the other hand, the resummation properties of the parton shower are retained. 
The difference in the actual Sudakov shape can entirely be attributed to the respective size and functional 
form of the phase-space boundaries. Further, within the region where the standard parton shower is able to 
fill the phase space (up to \t\ = /i|,), the L0(8)PS result scaled by a global factor K = (Tnlo/clo = 2.3 is in 
good agreement with the PoWHEG and Mc@Nlo with acut = 1 distributions. This is a consequence of the 
simple form of the virtual corrections resulting in a local ii'-factor B^^^/B, that largely coincides with the 
global X-factor. 

It was claimed in |25J that the difference between PoWHEG and Mc@NL0 predictions for Higgs-boson pro- 
duction originates from the local if-factor B/B in PoWHEG. Figure [2] shows explicitly that this claim is at 
odds with our findings. It displays predictions from a standard PoWHEG simulation and a modified one, 
where we replace B — >■ B. The same analysis is repeated with MC@Nlo at a = 1, except that in this case we 
replace B^^^ — B. It is manifest that the original and the modified PoWHEG results agree, except for a global 
X-factor, as was observed in [5S]. However, the same effect is found for MC@Nlo at a = 1. In contrast, 
both the modified PoWHEG and the modified MC@Nlo result differ in shape from the NLO prediction at 
high pt, while the Mc@Nlo result with a = 0.01 matches the NLO result exactly. The deviations between 
POWHEG and Mc@Nlo, which were observed in [25] can therefore not be attributed to the local if-factor 
B — > B. Instead they must originate from the unrestricted phase space in PoWHEG. 

Similar observations are made for the transverse momentum of the hardest jet of Fig [3] Here, the difference 
in the emission rate for H- and S-events, as discussed earlier, is plainly visible. Again, Mc@NL0 with 
occut = 0.01 agrees with the fixed-order NLO prediction for large jet transverse momenta. Also PoWHEG and 
MC@Nlo with Qfcut = 1 give similar results, overestimating the emission rate by more than 250%. 
This overestimation also feeds into the rapidity separation between the Higgs boson and the first jet in Fig.|4] 
The predictions of Mc@NL0 with oicut = 0.01 and Mc@Nlo with acut = 1 differ by a factor of two. The 
effects of decreasing acut are examined in Fig. 4a With decreasing acut, the rapidity difference distribution 
develops a pronounced dip, which becomes broader and flattens out until the fixed-order result is approached 
for acut ~0.01 and stabilises for even smaller acut- Again, this is explained by the different composition of 
the samples in terms of H- and §-cvcnts. Similar findings were reported in |39| . where tt-production was 
analysed and the Mc@NL0 method was compared to tree-level merging. The authors concluded that a dip 
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Rapidity separation betw. Higgs and leading jet in pp — )■ /r + X 




3 4 
Ay(h, 1st jet) 



(a) 



0.09 
0.08 
0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 



1-5 



0.5 



NLO 

MC@NLO a = 1 
MC@NLO a = 0.01 



Powheg 
LO®PS 
LOg)PS X 2.3 E 




3 4 
Ai/(h, 1st jet) 



(b) 



Figure 4: Rapidity separation between Higgs-boson and leading jet in inclusive Higgs-boson production 
{nih = 120 GeV) at E^^^ = 7 TeV. See Fig. [l] for details. 



in the leading jet rapidity spectrum, present only in the McQNlo prediction, was most likely an artifact 
of the incomplete phase-space coverage in Herwig's parton shower. While the dead zones of this parton 
shower are clearly different from the dead zones generated through ctcut in our approach, the resulting effect 
is surprisingly similar. It suggests that the dip, which is not present in either NLO fixed-order or tree- 
level merged results, must be attributed to exponentiation uncertainties. Fig. 4b shows a comparison with 
Powheg and LO(XiPS. Powheg predictions agree well with those from our Mc@NL0 with acut = 1 and 
from L0(8)PS with a global A'-factor, while MC@Nlo with acut = 0.01 agrees with standard L0(8)PS. 
Fig. [5] shows the rapidity spectrum of the Higgs boson, which can be defined at Born level and is thus 
described at NLO accuracy. No significant variation is observed in this spectrum when the parameter acut is 
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varied, as expected. One can thus conclude that both the Mc@NL0 and POWHEG techniques are consistently 
implemented in the event generator Sherpa. 

Summarising the above results, we find that the perturbative uncertainties associated with the PoWHEG 
method are large, due to the unrestricted phase space available for exponentiation. It is important to stress 
that this problem is in principle solved by the Mc@NL0 method and that we simply use our implementation 
of Mc@NL0 to quantify the efi'ect. 

4.2 Higgs-boson production in association with a jet 

This section presents predictions for the production of a Higgs boson in association with at least one jet. 
This process has not yet been investigated using either the Mc@NL0 or the PoWHEG approach. The Higgs 
mass is set to nih = 120 GeV. Jets at parton-level (i.e. before parton shower emissions take place) are defined 
using the inclusive k± algorithm [37] with R — 0.5 and p_L,min = 10 GeV. The independence of the results of 
the precise value of generation cut has been checked by varying p±.min = 5 ... 15 GeV. Again, the effective 
coupling of the Higgs to gluons, mediated by a top-quark loop, is modeled through an effective Lagrangian 
[35], and the virtual matrix elements here are computed by MCFM [ID]. 

Both the renormalisation and factorisation scales are set to the transverse momentum of the hardest parton- 
level jet in the event. The two powers of as in the effective gluon-gluon-Higgs coupling are evaluated at a 
fixed scale rrih, independent of the renormalisation scale /x. The reason for this choice is that scale variations 
of the effective coupling do not reflect a systematic uncertainty of the calculation at hand. They should 
be compensated by higher-order corrections to the effective coupling (i.e. diagrams with gluons attached to 
the top quarks in the loop), and will not be compensated by the one- loop corrections in the h+iet matrix 
elements that we use. Thus, while the respective scale uncertainty is part of the systematics of the final 
cross section prediction, its proper assessment must come from a different calculation. 

Our central prediction chooses Ofcut = 0.03 to minimise the uncontrolled exponentiation of non-logarithmic 
terms and distortions of the Sudakov shape of the Higgs boson's transverse momentum spectrum. The pre- 
dictions presented here include simulation of hadronisation |4T] and multiple parton interactions (MPI) [42] 
as well as hadron decays |^ and higher-order QED corrections to the h ^ tt decay [4^ . 

Uncertainty bands correspond to the variation of the renormalisation and factorisation scales in the range 
1/2/i . . . 2/z as well as to the variation of Ocut in a range from 1 to 0.001. Scales are varied only in the matrix 
elements, not in the parton shower. As before, only minimal cuts are applied in the analysis: two r-leptons 
with jryl < 3.5 andp_L > 25 GeV are required. Jets are defined using the inclusive A;_L-algorithm with R — 0.7 
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Figure 6: Higgs-boson transverse momentum for small values of p\ in inclusive Higgs-boson production 
(m/i = 120 GeV) at E^ms = 7 TeV. The contribution of S- and H-events in the MC@Nlo method 
is displayed separately for different values of acut- For acut = 1 and ckcut = 0-3 the sample 
consist primarily of S-events with small additions of negatively and positively valued H-events, 
respectively. 



16 



Transverse momentum of Higgs boson in pp ^ h -\- X 



MC@NLO a = 1 

H+S-events 

H-events only 

S-events only 




Sherpa 



30 



40 



50 



Transverse momentum of Higgs boson in pp ^ h -\- X 



^10- 









MC@NLO a = 0.3 








H+S-events _ 








H-events only = 




1 ^ 




S-events only 








■'"T--.^^; - - NLO 


~ Sherpa 

1 1 1 1 1 1 1 1 









30 



40 



50 



60 



70 



80 



90 100 

[GeV] 



Transverse momentum of Higgs boson in pp ^ h -\- X 



'U-i 




MC@NLO a = 0.03 






H+S-events _ 






H-events only : 






S-events only 






NLO 




"1-- "^t"'— 1_ 




~ Sherpa 













80 



10 20 30 40 50 60 70 

Transverse momentum of Higgs boson in pp — > + X 



90 100 
p'l [GeV] 



CI- 
TS 



-1- 


MC@NLO a = 0.003 




H+S-events _ 




H-events only - 




S-events only 




(absolute of S) 




NLO 


Sherpa I_ 





30 



40 



50 



60 



70 



80 90 100 

p'l [GeV] 



70 80 90 100 




p\ [GeV] 




Transverse momentum of Hig; 


5s boson in pp ^ h -\- X 




MC@NLO fl: = 0.1 


r ^T-L 


H+S-events _ 




H-events only = 




S-events only 




NLO 


~ Sherpa 





10 20 30 40 50 60 70 So 

Transverse momentum of Higgs boson in pp — > + X 



90 100 
[GeV] 





MC@NLO a - 0.01 




H+S-events _ 




H-events only : 




S-events only 




NLO 


Sherpa ! 

1 1 1 1 1 1 1 1 1 1 1 1 ' '1 1 1 





30 



40 



50 



60 



90 100 

p'l [GeV] 



Transverse momentum of Higgs boson in pp ^ h -\- X 









MC@NLO a = 0.001 




H+S-events _ 




H-events only - 




S-events only 




(absolute of S) 




"""•-L-u..^^^^ NLO 


Sherpa 





10 20 30 40 50 60 70 



90 100 

p'l [GeV] 



Figure 7: Higgs-boson transverse momentum for large values of p\ in inclusive Higgs-boson production 
(m/i = 120 GeV) at Ecms = 7 TeV. The contribution of S- and H-events in the MC@Nlo method 
is displayed separately for different values of acut- For acut = 1 and acut = 0-3 the sample 
consist primarily of S-events with small additions of negatively and positively valued H-events, 
respectively. 
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Figure 8: Prediction and uncertainties for the transverse momentum of the Higgs boson (left) and the 
leading jet (right) in Higgs boson plus jet events {rrih = 120 GeV) at i?cms — 7 TeV. 



and p_L_inin — 20 GeV. Their minimum transverse momentum requested in the analysis ensures that the 
parton-level event selection is inclusive enough to guarantee coverage of the full phase space. 

Figure |8] shows the transverse momentum spectra of the Higgs boson and the leading jet. Grey bands 
indicate the uncertainty due to the choice of ttcut , which will be referred to as exponentiation uncertainty in 
the following. Yellow bands show the scale uncertainty. It is manifest, that the exponentiation uncertainty 
is largest in the region of large transverse momenta, a fact that seems counter-intuitive at first. Despite the 
transverse momentum distributions being described at NLO, they can apparently be altered significantly by 



additional emissions in the parton shower. Taking a closer look at Eqs. (2.10l-(2.12), this fact becomes clear 



The proof of next-to-leading order accuracy of both the POWHEG and the Mc@NL0 formulae is based 
on an expansion of the Sudakov form factors to first order. Higher-order corrections do play an impor- 
tant role, however, especially in regions of the phase space where parton emission is logarithmically 
enhanced. This region is drastically enlarged if the upper integration limit of the exponentiated real 
emission is shifted from fip to s = E^j^g, as is the case in PoWHEG and Mc@NL0 with ctcut = 1 (cf. 
Sec. 4.1 1. This change alters the shape of the emission term in Eq. (2.10) compared to the emission 



term in Eq. (2.9 1. This finding is reinforced by the fact that in all distributions investigated, the 



predictions for acut = 0.03 coincide with those for Ucut — 0.001. 

The strong coupling is evaluated at a scale of the order of the relative transverse momentum between 
two partons in the emission term in Eq. (2.9), while it is evaluated at a different scale (here the 



transverse momentum with respect to the beam) in Eq. (2.11). This leads to an additional distortion 



of the spectrum of the real-emission term, which is formally a sub-leading effect, but which becomes 
important in the enlarged logarithmically enhanced regions of the emission phase space. 

The above effects are amplified in Fig. |9j as the observables there, the transverse momentum of the second 
jet and the scalar sum of all jets in the event, are described at leading-order accuracy only. 

In contrast, Figs. 10 and 11 show a rather mild dependence on acut in the shape of the distributions. 
Nonetheless, a change in the normalisation can be observed, which is explained by the fact that due to larger 
emission rates with increasing acut, the amount of radiative events is increased. Their kinematic distribution, 
however, does not seem to differ significantly. In fact, the large scale dependence observed in Figs. [8] and [9] 
is amplified here, especially for observables related to the hardest jet, indicating that, although /i — p± is an 
appropriate scale choice for this process, the canonical variation to ^ ~ leads to unphysical behaviour 
for jet rapidities beyond \r]\ « 3.4. 

As a consequence of the strongly varying jet rates with varying acut, the pseudorapidity separation between 
first and second hardest jet suffers from large uncertainties. It is displayed in Fig. [12] While the shape of 
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Figure 9: Predictions and uncertainties for transverse momentum of the second hardest jet (left) and the 
scalar sum of all jet transverse momenta (right) in Higgs boson plus jet events (mh = 120 GeV) 
at E,^, = 7 TeV. 




Figure 10: Prediction and uncertainties for the rapidity of the Higgs boson (left) and the leading jet (right) 
in Higgs boson plus jet events {mh — 120 GeV) at -Ecms — 7 TeV. 



the distribution is mostly unaffected, the normalisation varies strongly. The azimuthal separation between 
the two leading jets indicates how the exponentiation uncertainty affects the relative position of both jets in 
phase space. Compared to the pure next-to-leading order result, the back-to-back situation is amplified by 
a factor ^ 10 in the Mc@NL0 prediction. In such events, the Higgs-boson is likely to have been produced 
at p± 0, consistent with the observation of a depletion of events in the first bins of the Higgs transverse 
momentum spectrum for ofcut = 1, cf. Fig.|8] 



4.3 W-boson production in association with a jet 

In the following, W[—?' ej/]-f jet production in pp collisions at 7 TeV is studied. As for the case of Higgs-boson 
production, the Mc@Nlo implementation is validated against a fixed-order calculation. In addition the scale 
uncertainties are assessed by varying factorisation and renormalisation scales in Mc@NL0 by a factor of two 
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Figure 11: Predictions and uncertainties for the pseudorapidity (left) and angular (right) separation of the 
Higgs boson and the leading jet in Higgs boson plus jet events {ruh = 120 GeV) at Ecms = 7 TeV. 
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Figure 12: Predictions and uncertainties for the pseudorapidity (left) and azimuthal (right) separation of 
the two hardest jets in Higgs boson plus jet events {nih = 120 GeV) at i?cms = 7 TeV. 



in both directions. 

To make the validation as meaningful as possible effects from parton showering beyond the first emission in 
S-events, hadronisation and multiple parton scattering are not included in these parton level (PL) studies. 
Events are generated requiring at least OnG jet with p_\_,rain — 

10 GeV, defined according to the inclusive 
fcj_-algorithm [37 with R = 0.5. The transverse momentum of that hardest jet also sets the central value of 
factorisation and renormalisation scales. The level of exponentiation is fixed by acut = 0.03. 

The total generated cross sections of the NLO calculation, ctnlo = (4695±12) pb, and Mc@NL0, (Tmc@nlo = 
(4680 ± 30) pb, agree within statistical uncertainties. The total generated cross section in Mc@Nlo varies 
between a{pi/2) = (4950 ± 40) pb and a{fj. ■ 2) = (4510 ± 23) pb, i.e. up to 5.8% around the central value. 

Properties of the VF-boson, which is reconstructed from the truth neutrino and the lepton, are displayed in 
Fig. |13| Agreement between fixed-order and Mc@NL0 results in both the rapidity and transverse momentum 
spectra is found. Agreement is also found for the transverse momentum and pseudorapidity distributions of 



the electron, which are shown in Fig. 14 
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Figure 13: Transverse raomentum and rapidity of the Tl^-boson in > ev] + j production at the LHC. 
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Figure 14: Transverse momentum and pseudorapidity of the electron in W[-^ ev\ + j production at the 
LHC. 



In the event analysis jets are defined according to the inclusive fcj^ algorithm [37] with R — 0.7 and requiring 



p_L > 20 GeV. Properties of these jets are displayed in Fig. 15 The transverse momentum of the leading jet, 
an observable described at next-to-leading order accuracy, shows good agreement between both approaches. 
The scalar sum of jet transverse momenta {flr) is only described at leading-order accuracy and thus suffers 
from larger exponentiation uncertainties. This leads to a disagreement between the fixed-order result and 
Mc@NL0, up to the level of 20%. 

No distribution is affected significantly by scale variations beyond the change in total rate. This indicates 
that the functional form of the scale choice described at the beginning of this section works well for this 
process. 



5 Non-perturbative uncertainties and comparison to data 

In this section the Mc@NL0 method is further studied using VF-|-jet and 2'-t-jet production as a testing 
ground. Emphasis is now placed on the investigation of non-perturbative effects and associated systematic 
variations. 
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Figure 15: Transverse momentum of the leading jet and scalar sum of all jet transverse momenta in 
ev] + j production at the LHC. 




Figure 16: Transverse momentum and rapidity spectrum of the M^-boson in W[-^ ev] + j production at 
the LHC. 



5.1 Analysis of non-perturbative effects 

The Mc@NL0 method allows to generate fully hadronised events as an input for detector simulation or 
for direct comparison to measurements at the particle level. A question that naturally arises is, whether 
the theoretical uncertainties of the full Mc@NL0 simulation are then dominated by the perturbative or 
non-perturbative effects. We do not attempt to judge on this question here as it is not obvious on which 
grounds they can be compared at all, but rather point out that both non-perturbative corrections and 
non-perturbative uncertainties are modest compared to the intrinsic uncertainties of the parton-level result, 
which were investigated in the previous section. 

To this end, Mc@NL0 simulations for W-|-jet production are compared with a varying level of non-perturbative 
effects included: 

"Parton Level" 

Only the first emission off §-evcnts in McQNlo is generated in addition to the seed event. This is the 
same method that was used in the comparison to fixed-order results in Sec. |4.3| 

"Shower Level" 
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Figure 17: Transverse momentum and pseudorapidity spectrum of the electron in T4^[— > eiy] + j production 
at the LHC. 
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Figure 18: Transverse momenta and pseudorapidities of the two leading jets in W[-^ ev\ + j production 
at the LHC. 
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Figure 19: Correlations between PF-boson and leading jets in W[-^ eu] -\- j production at the LHC. 
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Figure 20: Hadronisation uncertainties for different observables studied in ei/] + j production at 

LHC. 
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All QCD emissions in the parton shower and QED emissions in the YFS approach are taken into 
account. 



"Shower+MPI" 

Multiple parton interactions and intrinsic transverse momentum of the beam hadron are additionally 
allowed for. 

"Hadron Level" 

Hadronisation and hadron decays are included to generate events at the particle level. 



All other event generation parameters have been chosen identical to Sec. |4.3[ 

Non-jet observables, like the rapidity and transverse momentum of the VF-boson and the rapidity and 
transverse momentum of the charged lepton are virtually unaffected by non-perturbative effects, as expected. 



This is shown in Figs. 16 and 17 



Jet observables are far more sensitive, as is exemplified by the transverse momentum and rapidity spectra of 
the leading and next-to-leading jet in Fig. [18] They show effects of "out-of-cone" radiation at shower level, 
which are partially compensated by the simulation of multiple scattering effects. Hadronisation again leads 
to softer jet spectra, such that the various effects can compensate each other. Although the tendency of the 
corrections stays the same, their precise magnitude depends on the jet algorithm and its parameters and 
will have to be investigated separately for each analysis. Figure [T9] exemplifies how non-perturbative effects 
can distort correlations between the two hardest jets. 

The uncertainties inherent in the hadronisation model were probed by switching from the Sherpa default 
cluster fragmentation model [41] to the Lund string model [45] in the implementation of Pythia [31, . Both 
models have been tuned to data from LEP and excellent agreement has been achieved. For this comparison 
to be meaningful the same perturbative events (i.e. identical random seeds) were subjected to the two 
hadronisation models. Effects from statistical fluctuations are thus cancelled. 



The differences were found to be negligible for the observables studied here, as displayed in Fig. 20 except 
for the specifically hadronisation-sensitive jet mass, where variations up to 20% occur. 



5.2 Z+jet production compared to Tevatron data 

Let us now turn to the comparison of Mc@NL0 predictions to data and a simple assessment of the systematic 
uncertainties associated with a significant part of the non-perturbative effects. Central predictions are made 
at the hadron level with Sherpa's default tune. An uncertainty band is generated by allowing for a variation 
of the MPI parameters within boundaries given by existing measurements of the underlying event. All event 
generation parameters have been chosen analogous to the ly-f jet case above. 

Figure [21] compares Z+jet production in the electron channel to a measurement from the CDF collabo- 
ration [46) . The reconstructed electrons are required to have transverse momenta p± > 25 GeV and an 
invariant mass of 66 < rUee < 116 GeV. Jets are defined using the midpoint cone algorithm [47] with R = 0.7 
and a split/merge fraction of 0.75. At least one jet with p± > 30 GeV and \y\ < 2.1 needs to be present and 
separated from both electrons by AR^j > 0.7. The cross section in the one-jet bin is predicted slightly too 
low. The cross section of the two-jet bin is significantly underestimated, as it is determined to leading order 
accuracy and subject to large uncertainties from the NLO-I-PS matching procedure as discussed in Sec. [4] 
The shape of the transverse momentum distributions agrees well with data. 

More characteristics of Z-boson plus jet production were investigated in a recent D0 analysis Events 
with two muons of invariant mass 65 < m^^ < 115 GeV and with at least one jet of p± > 20 GeV and 
\y\ < 2.8 were collected at a center-of-mass energy of 1.96 TeV. Jets were defined using the D0 midpoint 
cone algorithm [13] with R = 0.5 and a split/merge fraction of 0.5. Each jet had to be separated from both 



leptons by AR^j > 0.5. A comparison of Monte-Carlo predictions with this measurement is shown in Fig. 22 
The agreement is fair, except for the pj^-spectrum of the first jet, where a deficiency of the Monte-Carlo 
result at intermediate and high p± is observed. 

A further measurement of Z+jet production in the electron channel was presented by the D0 collaboration 
in [50]. Each electron is required to have p± > 25 GeV and the mass window 65 < niee < 115 GeV is 
enforced. Jets are defined using the midpoint cone algorithm [49 with R = 0.5 and a split/merge fraction 
of 0.5. At least one jet with p± > 20 GeV and \ri\ < 2.5 must be present in the event. Experimental data 
were normalised to the inclusive Drell-Yan cross section. This quantity is not predicted by the Mc@NL0 
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Figure 21: Inclusive jet cross sections and inclusive transverse momentum distributions of all and all-but- 
the-hardest jets compared to CDF data |46) . 



simulation of Z+jet, and therefore the Mc@NL0 results are scaled by a global factor of 0.35 such that 
the normalisation of the one-jet-rate agrees with data. Figure 23 displays the comparison of these scaled 
MC@Nlo predictions with the results from D0 . The transverse momentum shape of the leading jet is 
well described by the simulation. The two-jet rate, described here at leading-order accuracy, seems to be 
underestimated compared to data, but the shape of the sub-leading jet's p± spectrum is relatively flat. For 
the third jet, which is generated in the parton-shower approximation, both the rate and the shape of the 
spectrum are not simulated correctly. 

To quantify the success of the next-to-leading order calculation it is important not only to investigate single- 
particle spectra, but also correlations between the Z-boson and the hardest jet. They might give some insight 
into the genuine one-loop effects in Z+jet production. Therefore, the analysis strategy of a measurement 
presented by the D0 collaboration in [51^ is pursued. Opposite-sign muons with p± > 15 GeV and an 
invariant mass of 65 < m^^; < 115 GeV are required in association with at least one jet of p± > 20 GeV and 
\y\ < 2.8. Jets are defined using the midpoint cone algorithm [15] with R = 0.5 and a split/merge fraction of 
0.5. Two event samples are defined, one requiring the transverse momentum of the reconstructed Z boson 
to be above 25 GeV, the other requiring it to be above 45 GeV. Figure [24| compares MC@Nlo predictions 
to the measurement. The shape of the rapidity distributions is matched fairly well, again with the total 
rate prediction at the lower end of the uncertainty band. For the azimuthal correlation distribution the 
shape shows significant deviations from data: Only the back-to-back configuration is described well, but for 
A(j) < IT the Mc@NL0 prediction underestimates the data. This signals the uncertainty related to emissions 
beyond the hardest one, which are only generated at leading-order or parton-shower accuracy and which are 
also subject to large NLO-I-PS matching systematics as discussed in Sec. |4] 



5.3 W+jet production compared to LHC data 

The production of a boson in association with at least one hard jet in proton-proton collisions at -^s = 7 
TeV has been studied by the ATLAS cohaboration at the CERN LHC [lU] . 

In the electron channel events are selected by requiring an electron with p± > 20 GeV defined at the 
particle level to include all photon radiation within a Ai? = 0.1 cone. Only electrons in the fiducial volume 
\r]\ < 1.37 or 1.52 < \r]\ < 2.47 are taken into account. E™^^^ at the particle level has been defined through 
the leading neutrino in the event which is required to have p±^i, > 25 GeV. The transverse mass cut is placed 

at TOT = ^J'^pip'i (1 - cos{<j)e -(f>^))> AO GeV. 

Jets are reconstructed using the anti-fcf algorithm with R = 0.4 and have been taken into account if p'^ > 
20 GeV, |?7jet| < 2.8 and Ai?(^,jet) > 0.5. Muons, neutrinos and the leading electron were excluded from 
the input of jet reconstruction. 

The comparison in Figure [25] shows good agreement of the Sherpa hadron level prediction with ATLAS 
data in the shapes of differential distributions on the one hand, and discrepancies in the prediction of the 
total rate on the other hand. While the one-jet rate is still predicted fairly well, the two and three-jet rates 
are significantly too low. The latter are only predicted at leading-order and parton-shower accuracy and 
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Figure 22: Transverse momentum and rapidity distributions of the reconstructed Z boson (top row), and 
the total cross section, the transverse momentum and rapidity distributions of the hardest jet 
(bottom row) in Drell-Yan production in association with at least one jet compared to D0 data 



also subject to the inherent uncertainties from NLO+PS matching. 



6 Conclusions 

In this publication, the MC@Nlo and POWHEG methods to match NLO QCD matrix elements with parton 
showers have been compared, with a special emphasis on some issues that generate large differences between 
their predictions. In particular, these issues are related to 

• the treatment of sub-leading colour configurations in Mc@NL0 

• the choice of scales in PoWHEG and Mc@NL0 and 

• the exponentiation of non-leading terms in PoWHEG. 

Before discussing the findings related to these issues, it is worth summarising some other results presented 
in this publication. They refer to the impact of scale variations at the parton-level of the processes studied 
here, and to the effect of hadronisation and the underlying event on a number of observables: 
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Figure 23: Transverse momentum distributions of the three hardest jets in Drell-Yan production in asso- 
ciation with at least one jet compared to D0 data [SO]. All SherpA predictions are scaled with 
a common factor to account for the unknown normalisation to the inclusive process. 
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Figure 24: Azimuthal and rapidity difference distributions of the reconstructed Z boson and the hardest jet 
and their rapidity average in Drell-Yan production in association with at least one jet compared 
to D0 data [5T] . 



• Scale uncertainties 

In W/Z{+iets) production, uncertainties due to scale variations are typically about 5-10%, but they 
can increase up to about 25% in the high-p_L region of jet production, where higher jet multiplicities 
become relevant. In Higgs production processes scale uncertainties tend to be significantly larger, by 
at least a factor of 2, which is indicative of large higher-order corrections. 

• Underlying event corrections and uncertainties 

Multiple parton interactions tend to increase jet production rates, and they have a particularly large 
influence on observables which are sensitive to high jet multiplicities, like Ht- They typically do not 
affect inclusive observables, with changes of about 5% or below. Their associated uncertainty is small. 

• Hadronisation corrections and uncertainties 
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Figure 25: Inclusive jet cross sections and transverse momentum spectra of the hardest and second hardest 
jet in ei/] + j production at the LHC compared to ATLAS data [TU] . 



Hadronisation corrections are similar in magnitude to parton-shower and multiple-scattering effects. 
The associated uncertainties are about 5%, except for observables such as jet masses, where they can 
increase to about 20%. 

Let us now turn to the uncertainties induced by the NLO+PS matching. 

Regarding the problem of infrared divergences in sub-leading colour configurations, it should be noted that 
up to now, only one critical process - heavy quark pair production - has been implemented in a publicly 
available and fully documented event generator. In this case the problem was solved by introducing a factor 



that modifies the integral in Eq. (2.111 such that contributions from the region of soft-gluon emission are 
removed. A similar strategy seems to be employed in the recently presented aMc@NLO code, but it has 
not been published or made publicly available yet. It is therefore not yet clear whether this technique is 
sufficiently process-independent to allow for a general implementation of Mc@NL0. In this publication a 
different approach has been followed, which aims instead at constructing an exact solution, irrespective of 
the process considered. The basic idea is to correct the Casimir operators used in the parton shower such 
that the full colour structure in the soft-gluon limit is obtained. This introduces both positive and negative 
splitting functions, which is unproblematic for NLO calculation but implies "anti-probabilistic" evolution in 
a parton shower. To resolve this issue, a weighting procedure has been introduced. It is anticipated that this 
modification will eventually allow the systematic inclusion of sub-leading colour terms in the parton shower. 

The exponentiation of non-leading terms in POWHEG is an issue that has been known for some time. Actual 
differences in how PoWHEG and Mc@NL0 predict the p± spectrum of the Higgs boson are documented 
in the literature |25j . They have typically been attributed to a different treatment of higher order, in 
particular NNLO, corrections (see Sec. 4.3 of [25]). We find that the difi^erence is related to exponentiation 
uncertainties, naturally affecting the next-to-next-to leading order, but continuing to all orders. The upper 
scale in the logarithms resummed in PoWHEG is related to the total hadronic centre-of-mass energy, which 
is substantially different from the scale typically used in the resummation programme |521 153j , where a value 
of the order of the Higgs boson mass is assumed. It is thus clear that a suitable constraint on the emission 
phase space in POWHEG must be defined. This problem has already been solved by the Mc@NL0 method. 
On a related note, the occurrence of the well-known dip around zero in the Mc@NL0 simulation of the 
Vj ~ Vh distribution, the rapidity difference of Higgs boson and accompanying jet, has been studied here. 
Cutting on the phase space of the Catani-Seymour subtraction kernels and refiecting these cuts in the parton 
shower, as enforced by Mc@NL0, yields dead zones of parton emission in the latter. While these dead zones 
are clearly different from the dead zones in the HerwiG parton shower, the resulting effect is surprisingly 
similar. In fact, both including the full phase space and comparably tight cuts on the phase space leads to 
a vanishing dip; this finding suggests that the dip, which is not present in the NLO calculation, must be 
attributed to true exponentiation uncertainties. 

To complete our implementation of the Mc@NL0 method in the sense that only logarithms related to 
emissions on scales below the factorisation scale are exponentiated, physically more meaningful cuts on the 
phase space of Catani-Seymour subtraction kernels are currently investigated. Devising a similar technology 
for PoWHEG would result in two formally completely equivalent algorithms, but imply a larger computational 
effort for event generation in PoWHEG due to the additional numerical integration in B and the matrix- 
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element correction of the weighted parton shower, see [T3]. For this reason, in favour of the Mc@NL0 
method, the further development of POWHEG methods in the Sherpa framework will be abandoned. 
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A Comment on the NLL accuracy of the PowHEG formula 

In [ST it was proven that the POWHEG formula yields NLL accurate predictions, if the strong coupling is 
modified according to |55| . The proof rests on the crucial assumption that when logarithms of the form 
L = log(pj^/(5) are resummed, the dependence on the resummation scale Q introduces corrections beyond 
next-to-leading logarithmic order. This argument holds as long as Q differs from the hard scale in the 
process, say mt in Higgs-boson production, by a factor of order one [SHI [511 IMl [S7j ■ As pointed out in [S^ 
large spurious subleading logarithms may otherwise be introduced, that are associated with the overall 
normalization of observables, rather than radiative corrections. In other words, in the large-fcy region the 
resummed results lose their predictivity and should be replaced by conventional fixed-order calculations [57j . 
Following these arguments, the reasoning in ,54, is only correct as long as the resummation scale and the 
hard scale are sufficiently similar. This is not the case in the PoWHEG method when it is applied to hadronic 
collisions where Q — >■ Ecms, the total hadronic centre-of-mass energy, while the hard scale is of order rrih, for 
example. 
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